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The Saffman- Taylor viscous fingering problem is investigated for the displacement of a non- 
Newtonian fluid by a Newtonian one in a radial Hele-Shaw cell. We execute a mode-coupling 
approach to the problem and examine the morphology of the fluid-fluid interface in the weak shear 
limit. A differential equation describing the early nonlinear evolution of the interface modes is de- 
rived in detail. Owing to vorticity arising from our modified Darcy's law, we introduce a vector 
potential for the velocity in contrast to the conventional scalar potential. Our analytical results ad- 
dress how mode-coupling dynamics relates to tip-splitting and side branching in both shear thinning 
and shear thickening cases. The development of non-Newtonian interfacial patterns in rectangular 
Hele-Shaw cells is also analyzed. 

PACS numbers: 47.50.+d, 47.20.Ma, 47.54.+r, 68.05.-n 

I. INTRODUCTION 

The celebrated Saffman- Taylor instability Q arises when a less viscous fluid pushes a more viscous one in the thin 
gap of a Hele-Shaw cell. The less viscous fluid can be either injected at an end of a channel-shaped cell (rectangular 
geometry) , or from the center of the cell (radial geometry) [|j . In both geometries the interface separating the fluids 
may deform, leading to the formation of hngerlike patterns, commonly known as viscous fingers. During the last four 
decades, the viscous fingering instability has been extensively studied, both theoretically Q and experimentally [Q. 
Much of the research in this area has examined the case in which the fluids involved are Newtonian. For Newtonian 
fluids, it is observed that the fingers grow and compete dynamically, resulting in a single stable finger in the rectangular 
geometry, and in patterned structures markedly characterized by the spreading, and subsequent splitting of the finger 
tips in the radial setup. 

A whole different class of interfacial patterns arise when the Saffman- Taylor instability is studied by taking the 
displaced fluid as non-Newtonian [Q, |^| . In contrast to most Newtonian fluids, non-Newtonian fluids differ widely in 
their hydrodynamic properties, with different fluids exhibiting a range of effects from elasticity and plasticity to shear 
thinning and shear thickening. Experiments using non-Newtonian fluids in radial and rectangular Hele-Shaw cells have 
revealed a wide variety of new patterns, showing snowflake-like shapes ^, || and fracturelike structures [Q, |^] . Instead 
of the traditional, tip-splitting-dominated Newtonian patterns, these experiments |^, 0, ||] exhibit dendritic fingers 
and side branching. This morphological diversity and rich dynamical behavior motivated a number of theoretical 
studies of the problem §, |o|, [H], [if, ||. One major difficulty faced by researchers is that, as opposed to the 
Newtonian case, the pressure field is no longer Laplacian. This implies a serious theoretical challenge, since a priori 
one would not be allowed to directly apply a Darcy's law approach to attack the problem. 

The instability of radial Hele-Shaw flows involving non-Newtonian fluids has been studied theoretically by Sader 
et al. ||. They considered power law fluids, and performed a linear stability analysis without making use of Darcy's 
law. Essentially, they showed that decreasing the power law index dramatically increases the growth rates leading 
to a more rapid development of the fingering patterns. A Darcy's law- type approach has been proposed by Bonn 
and collaborators JlO| . Bonn et al. suggested a modified Darcy's law including a shear rate dependent viscosity, and 
showed that, within their approach, the pressure field remains Laplacian. In a series of interesting papers Kondic et 
al. [[n} [lj] , and subsequently Fast et al. |1^] extended reference [[to) ideas and derived a generalized Darcy's law from 
first principles, where viscosity depends upon the squared pressure gradient. It turns out that the Darcy's law formula 
proposed by Bonn et al. jl0| follows from the more basic version rigorously derived by Kondic et al. [[llj [lj] . Efficient 
numerical simulations performed in references ]l2| , |l3| have shown that shear thinning can suppress tip-splitting and 
leads to the formation of dendritic structures, presenting a clear side branching behavior. 
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Theoretical studies of the fully nonlinear stages of Hele-Shaw flow with non-Newtonian fluids rely heavily on 
intensive numerical simulations jl2|, |l3). On the analytical side, the structure of the fingering dynamics in such 
complex fluids is largely restricted to linear stability investigations ||, [ll| . Much less attention has been paid to the 
analytical investigation of the dynamics that bridges the initial (purely linear) and final (fully nonlinear) time regimes. 
In addition, theoretical as well as experimental analyses of flow of shear thickening fluids in Hele-Shaw cells still need 
to be addressed. In this paper we carry out the analytical weakly nonlinear analysis for the intermediate stages 
of evolution, and examine both shear-thinning and shear thickening cases. We adapt a weakly nonlinear approach 
originally developed to study Newtonian Hele-Shaw flows |l4|, |lq] , to the non-Newtonian situation. We focus on 
the onset of the nonlinear effects, and try to understand how mode-coupling dynamics leads to basic morphological 
features and behaviors observed in non-Newtonian Hele-Shaw flows. 

The article is organized as follows: section [n] formulates our theoretical approach. We perform a Fourier decom- 
position of the interface shape, and from an alternative form of Darcy's law study the influence of weak shear effects 
on the development of interfacial patterns. In contrast to the analysis for Newtonian fluids, conventionally based 
on a scalar velocity potential, we employ a vector velocity potential capable of describing vorticity arising from the 
non-Newtonian fluid flow. Coupled, nonlinear, ordinary differential equations governing the time evolution of Fourier 
amplitudes are derived in detail. Section IH discusses both linear and weakly nonlinear dynamics. It conce ntrates 
on the effect of shear thinning and shear thickening on finger tip-splitting and side branching. Section III A briefly 
discusses our linear stability results. Line ar res ults are useful and instructive, but do not allow accurate predictions of 
important interfacial features. In section III B we show that some of such features can indeed be predicted and more 
quantitatively explained by our analytical mode-coupling approach. At second order we find a mechanism for finger 
tip-splitting in non-Newtonian Hele-Shaw flow: it is suppressed (favored) for shear thinning (thickening) fluids. Our 
results indicate absence of side branching in the weak shear limit a nd ea rly flow stages, but suggest that it could be 
enhanced (inhibited) for shear thinning (thickening) fluids. Section III C| discusses mode-coupling in rectangular flow 



geometry. Our chief conclusions are summarized in Sec. [V 




non— Newton ian 
fluid 



FIG. 1: Schematic configuration of radial flow in a Hele-Shaw cell. The inner fluid is Newtonian and has negligible viscosity. 
The outer fluid is non-Newtonian. All physical parameters are defined in the text. 



II. THE MODE-COUPLING DIFFERENTIAL EQUATION 



The Hele-Shaw cell (Fig.l) consists of two parallel plates separated by a small distance. The cell thickness b is 
considered to be much smaller that any other length scale in the problem, so that the system is essentially two- 
dimensional. Consider the displacement of a viscous, non-Newtonian fluid by a Newtonian fluid of negligible viscosity 
in such confined geometry. The surface tension between the fluids is denoted by a. The Newtonian fluid is injected 
at a constant areal flow rate Q at the center of the cell. 

The initial circular fluid-fluid interface is slightly perturbed, 1Z = R(t)+( (9, t) ((,/ R <C 1), where the time dependent 
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unperturbed radius is given by 



R(t) = JB* + 9±, 



(1) 



Rq being the unperturbed radius at t — 0. The interface perturbation is written in the form of a Fourier expansion 



C(0,*)= J2 C«(*)exp(m0), 



(2) 



where Cn(*) = (1 /27r) £(9,t) exp (—in9) d9 denotes the complex Fourier mode amplitudes and n—0, ±1, ±2, ... 
is the discrete polar wave number. In our Fourier expansion (|^) we include the n = mode to keep the area of the 
perturbed shape independent of the perturbation £. Mass conservation imposes that the zeroth mode is written in 
terms of the other modes as Co = — {1/2R) \Cn(t)\ 2 . 

The relevant hydrodynamic equation for Newtonian Hele-Shaw flows, Darcy's law |l], ||, states that the gradient of 
the pressure is proportional to the fluid velocity, and oriented in the opposite direction with respect to the fluid flow 



Wp - 



12(i 



(3) 



where v = v(r, 9) and p = p(r, 9) are, respectively, the gap-averaged velocity and pressure in the fluid. The viscosity 
of the fluid is represented by fx. To model non-Newtonian Hele-Shaw flows, we use a suitable form of the Darcy's law 
in the weak shear thinning or thickening limit. We follow Bonn et al. |To[ and consider their shear rate dependent 
viscosity in the weak non-Newtonian limit (wt) 2 <C 1 



H(w 2 ) » /i [l - (1 - ol) u> 



VI 



(4) 



where lo ~ v/b is the shear rate, v = |v|, r denotes the characteristic relaxation time of the fluid, and fio is a constant, 
zero-shear viscosity. The parameter a measures the shear dependence: a = 1 corresponds to the Newtonian fluids, 
and a < 1 (a > 1) gives the shear thinning (thickening) case. 

By substituing Eq. (g) into Eq. (g) we obtain an alternative form of Darcy's law ideally suited to describe weak 
non-Newtonian effects 



VP = -v 



5v 2 w. 



(5) 



Here P = [& 2 / (12/j,o)] P represents a generalized pressure field, and 5 = (1 — a)(r/6) 2 is a small parameter that 
expresses the non-Newtonian nature of the displaced fluid: S — corresponds to the Newtonian case, while 5 > 
(5 < 0) describes the shear thinning (thickening) case. Taking the divergence of Eq. (||) and using incompressibility 
(V • v = 0) the pressure is seen to be anharmonic (nonvanishing Laplacian). Our subsequent analysis incorporates 
this effect, in contrast to the treatment in Ref. |Iofl . Indeed, an alternate route to our Eq. (|5|) is to start with viscosity 
depending on the square of the pressure gradient (as in Kondic, et al. fTT| ), then approximate the pressure gradient 
with the velocity for small S. 

Our pertubative approach keeps terms up to the second order in £ and up to first order in S. Considering the fact 
that the viscosity of the Newtonian fluid is negligible, the generalized pressure jump condition at the interface can be 
written as [HI 



P\n = -J k 



II l TC > 



(6) 



where 7 = b 2 a / '(12 hq) and ac || is the curvature in the direction parallel to the plates. 

The weakly nonlinear approach to radial, Newtonian Hele-Shaw flow developed in reference [15), related the fluid 
velocity to a scalar velocity potential v = — V0, this replacement made possible by the irrotational nature of the flow 
for Newtonian fluids. For non-Newtonian fluids, in contrast, flows governed by the modified Darcy's law (JsJ) exhibit 
vorticity. Hence we perform our calculations using a vector potential v=VxA. The most general form of the vector 
potential can be written as 



2n ^— ' \ r 



exp [inb 



(7) 
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where A mn are the Fourier coefficients of the velocity vector potential and z is the outward unit-normal to the upper 
cell plate. The radial and polar components of the fluids velocities are 



Q_ 

2wr 



inA ri 



R" 



j I exp (mb 



and 



ug = mA mn ( —^T[ ) exp (in9). 



(8) 



(9) 



Note that the vector potential reduces to the unperturbed steady flow with a circular interface (v r — Q/2nr, vg = 0) 
when R — > and also when r — > oo. 

We exploit the fact that VP must be curl free, and impose the so-called solvability condition V x VP = 0. It 
simplifies the general form of the vector potential expansion given in Eq. (Q). The solvability condition reveals that, 
without loss of generality, one can rewrite the vector potential as 



^ njtO V 7 



exp (inQ) + S 



/ J?\ 1 | 

n/0 V 7 J J 



(10) 



replacing the array of coefficients A mn with the simpler set of A n and B n . Observe that the vector potential (|T^) is 
simply a superposition of a purely Newtonian term (oc 5°, coefficients A n ) and a non-Newtonian contribution (cx 6 , 
coefficients B n ) 



A = A 



(11) 



The flow described by Ajv is irrotational, while AjvAf has a curl. 

Similarly, we express the pressure of the outer fluid as a sum of Newtonian and non-Newtonian pressures, and 
propose a general form for their Fourier expansion 



P = P N + P NN , 



where 



P 



N 



Q , (T 



n^O 



R 



\n\ 



exp (in6), 



(12) 



(13) 



and 



P 



NN 



QV 1 , W 1 

7 n^O V 7 



(in0) 



(14) 



The gradient of the complex pressure field (jT^) must satisfy the non-Newtonian Darcy's law given by Eq. (|5|). By 
inspecting the r and 9 components of and by examining the Newtonian and non-Newtonian components of it, we 
can express the Fourier coefficients of P^r, Pnn and K^n in terms of the Fourier coefficients of A^r, 



p n = sgn(n) (iA n ) , 



(15) 



q n = -iA n (3{ri) 



n ( Q 



Q 



\n\ \2tt J 4tt (|tt.| + 1) 



— ^2 mA m (n - to) A n - m k (n, to) , 



(16) 



B n = -A n a(n) I — 



Q 



2tt/ 4tt (ln.1 + 1) f— 1 , 



to (iv4 m ) (n - to) A„_ m /i (n, to) , 



(17) 
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where in order to keep the results in a compact form, we introduced the coefficients 

H(H-i) 



a (n) 



2(|n| + 1) 



(18) 



(3{n) 



H(|n|+3) 
2(|n| + l) ' 



(19) 



and 



h (n, to) = (|n| + 2) sgn(n — to) — — {3 — sgn [to (n — to)]} . 



fc (n, to) = — rt sgn(n — to) + — (\n\ + 2) {3 — sgn [to (n — ni)]} . 



(20) 



(21) 



Note that sgn(n) = 1 if n > and sgn(n) = — 1 if n < 0. 

Using Eqs. (|ll|)-([l7|), which are consistent with the solvability condition and Darcy's law (||), we can derive the 
general expression of the vector potential Fourier coefficients in terms of the perturbation amplitudes. To fulfill this 
goal, consider the kinematic boundary condition p, H 



an 



ion . . 

r o0 



(22) 



which states that the normal components of each fluid's velocity at the interface equals the velocity of the interface 
itself. By expanding Eq. ( p2] ) up to the second order in £ and up to first order in 8 we find the coefficient of the vector 

potential corresponding to the n-th evolution mode, An \ in terms of 5 and the fc-th order in £ (fc = 1, 2) 



R- R 

Sn i Sri 

n n 



1 + 5R 2 a (n) 



SRR 



h OR u(n, to) 



2(|n| + 1) 



1 |to| .« . 

— I h di? u(n, to) 

n to 



2(|n| + 1) 



— ^ h(n,m) CmCn-r 



(23) 



where the overdot denotes total time derivative, R = Q/(2irR) represents the unperturbed interface velocity, and the 
coefficients 



, . \m\ a(m) h(n,m) 
u (n, to) = a(n) J-i - 2^^ - [ , 
to ml (\n\ + 1) 



(24) 



, 1 |to|\ afm) h(n.m) 
v(n 7 m)=a(n) - + O - 

,n to / to 2 (|n| + 1) 



(25) 



To conclude our derivation we need one more step. The vector potential coefficients can be introduced into the 
pressure jump condition (Eq.(^)) and using Darcy's law (Eq.(^)) one can finally find the equation of motion for 
perturbation amplitudes £„. We present the evolution of the perturbation amplitudes in terms of 6 and the fc-th order 
in the perturbation amplitude £ 



(26) 



where 



C« = A(n)Cn, 



(27) 
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A(n) = |(|n|-l)--^|n| 



(n 2 - 1) + SR 2 



R 



R 



\n\-l 2 7 \n\ in 2 - l)~ 
Inl + 1 _ uj + j 



(28) 



is the non-Newtonian linear growth rate first reported in Kondic, ei aZ. |0] , and 



Ci 2) = X! [ f JvK m ) + <5 f mK m )] CmCn-m+ [Gn (n, to) + S G NN (n, to)] CmCn-m 

+ <5 H NN(n,m) CmCn-m + 5 X! J NN(n,m) ( m (n-m- (29) 

In Eq. ( p9| ) the coefficients Fn,Fnn, Gn, Gnn, Hnn, and Jatat represent the second order Newtonian (N) and non- 
Newtonian (NN) terms. Their detailed functional form are presented in the appendix. An important feature of the 
second order coefficients is that they present special reflection symmetries 

C (n, —to) = C {—n, to) 

C (— n, — to) = C (n, to) (30) 

where C — Fjf, fjvjvj ^JV; Gnn, Hjy n , and Jatat, respectively. Equation ( |2^ ) is the mode-coupling equation of the non- 
Newtonian Saffman- Taylor problem in radial geometry. It gives us the time evolution of the perturbation amplitudes 
Cn; accurate to second order, in the weak shear limit. The rest of the paper uses Eq. ( p6| ) to study the development 
of interfacial instabilities, and to examine how the non-Newtonian parameter 6 affects pattern morphology. 



III. DISCUSSION 



In the next three subsections we use our mode-coupling approach to investigate the interface evolution at first and 
second order. To simplify our discussion it is convenient to rewrite the net perturbation (|^) in terms of cosine and 
sine modes 

oo 

C(0, t) = Co + M*) cos (^) + b n{t) sin(nfl)] , (31) 

n=l 

where a n — ( n + £-n and b n = i (£„ — C-n) ar e real-valued. Without loss of generality we may choose the phase of 
the fundamental mode so that a n > and b n = 0. 



A. First order 



Although at the level of linear analysis we do not expect to detect or rigorously predict important nonlinear effects 
such as tip-splitting and side branching, some useful information may still be extracted. A nice example of how 
purely linear results can help to understand complicated morphological features appearing in non-Newtonian radial 
Hele-Shaw flows is found in references JO], |l2). By studying their linear growth rate, Kondic et al. [[y| [lj] found that 
shear thinning decreases the wave number of maximal growth, increases the maximum growth rate, and tightens the 
band of unstable modes. Based on this increased selectivity of wavelengths, they postulated that shear thinning can 
lead to suppression of tip-splitting. Their speculations have not been further investigated analytically, but instead 
have been supported by their own intensive numerical simulations Jl2|, |lF 



28J), agrees with Kondic et al. |llj equivalent 
to gain insight into shear thinning/thickening 



Our first order results, namely our linear growth rate expression 
formula, if we set their constant C equal to R 2 . We use growth rate (E 
behavior. Figure 2 plots A(n) as a function of mode number n for three different values of the non-Newtoniam 
parameter: a) S = 0; b) 5 = 0.05; and c) 8 = —0.05. By inspecting Fig. 2, we notice that, unlike the shear thinning 
case (6 > 0) discussed in references |ll| |lj] , shear thickening (6 < 0) widens the band of unstable modes, and decreases 
growth rate for the wave number of maximal growth. So, by applying similar arguments as those used by Kondic et 
al. [p"H [l2f for the weak shear thinning case, we postulate that shear thickening would lead to enhanced tip-splitting. 
Unfortunately, for the shear thickening radial Hele-Shaw flow, both numerical simulations and experiments are not 
available in the literature to confirm this claim. Of course, it is also of interest to study such a possibility analytically. 
This is one of the topics we examine in Sec. IIIB. 
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FIG. 2: Linear growth rate ( |28[ ) as a function of mode number n for three different values of the non-Newtoniam parameter: 
5 = (black solid curve), +0.05 (black dashed curve), -0.05 (gray dashed curve). The injection rate Q = 2n cm 2 /s, R = 1 cm, 
and 7 — 1/200 cm 3 /s. The units of A(n) and 5 are s _1 and (cm/s) -2 , respectively. 



B. Second order 



On the analytical side, one must go beyond linear analysis in order to investigate in more detail the main morpho- 
logical features found in Hele-Shaw flows with both shear thinning and shear thickening fluids. To this day, besides a 
few first order linear growth studies Ell |l2|, |l3) , there are no other analytical results concerning non- Newtonian fluid 
evolution in radial Hele-Shaw geometry. If on one hand the non-Newtonian problem is vastly richer, it is also true 
it is much less amenable to analytic attack than its Newtonian counterpart. Our mode-coupling approach intends to 
provide a better and more complete analytical understanding of the complex non-Newtonian time evolution dynamics 
in Hele-Shaw cells. 



1. Action of the parameter 3 on finger tip-splitting 

We use the mode-coupling equation ( |26| ) to investigate the coupling of a small number of modes. At second order 
the most noteworthy effect refers to the action of the non-Newtonian parameter 5 on finge r tip -splitting. Tip-splitting 
is related to the influence of a fundamental mode n on the growth of its harmonic 2n For consistent second 

order expressions, we replace the time derivative terms a n and b n by A(n) a n and A(n) 6„, respectively. Under these 
circumstances the equation of motion for the cosine mode 2n is 

a 2n = A(2n) a 2n + ^ T(2n, n) a? n , (32) 

where the function that multiplies a^, T(2n, n), is called the tip-splitting function and its general expression is 

T(n,m) = FN(n,m) + \(m) GN(n,m) + 6 [FNN(n,m) + \(m) GNN(n,m) 

+ \(n — m) HNN(n,m) + \(m) \(n ~ m) JNN(n,m)]. (33) 

Equation ([52]) shows that the presence of the fundamental mode n forces growth of the harmonic mode 2n. The 
function T(2n, n) acts like a driving force and its sign dictates if finger tip-splitting is favored or not by the dynamics. 
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If T(2n,n) < 0, ai n is driven negative, precisely the sign that leads to finger tip-widening and finger tip-splitting. If 
T(2n, n) > growth of a% n > would be favored, leading to outwards-pointing finger tip-narrowing. 




FIG. 3: Variation of T(2n, n) as a function of the non-Newtonian parameter 5, for three different Fourier modes (n = 4, 5, 6). 
The injection rate Q = 2-k cm 2 /s and 7 = f/IOO cm 3 /s. The units of T(2n, n) and 5 are (cm s) _1 and (cm/s) -2 , respectively. 



To investigate the influence of the non-Newtonian parameter S on tip-splitting behavior at second order, we plot in 
Fig. 3 the behavior of T(2n, n) as a function of 6, for a few Fourier modes (n — 4, 5, 6). To simplify our analysis we 
adopt an instantaneous approach: we consider a particular R and R combination, using the identity RR = Q/2ir, at 
the onset of growth of mode 2n [using the condition A(2n) = 0] in the Newtonian limit 5 = 0, where we know T(2n, n) 
is negative |T^]. We see from Fig. 3 that the curve for mode n = 4 lies below the curves associated to n = 5, 6. This 
is an expected behavior since smaller values of n would mean more room for the existing fingers to split. We also 
observe that the curves associated to smaller n show a stronger angular inclination with respect to the horizontal 
(5-axis. Therefore, lower Fourier modes would be more sensitive to variations in S. 

Further inspection of Fig. 3 reveals that, for 5 > 0, T(2n, n) becomes less negative as 5 increases, meaning that 
the interface has less tendency toward tip-splitting. In contrast, for 8 < we observe that T(2n,n) becomes more 
negative as the magnitude of 6 < increases, indicating an enhanced tendency of the fingers to split at their tips. We 
conclude that tip-splitting is suppressed for shear thinning fluids, and enhanced for shear thickening ones. This fact 
is more clearly illustrated in Fig. 4. In Fig. 4 we plot the fluid-fluid interface for a certain time (t = 30 s), considering 
the interaction of two cosine modes (a fundamental n = 4 and its harmonic 2n = 8), for three diferrent values of 5: 
(a) 6 = (solid curve); (b) 6 > (black dashed curve); and (c) S < (gray dashed curve). From Fig. 4 it is evident 
that there is a stronger splitting in the shear thickening case. 



These second order effects regarding tip-splitting are consistent with the first order effects described in section [II A . 
With respect to shear thinning behavior, our analytical results agree with numerical simulations fL2"| , [f3| and exper- 
iments [M o] of fully nonlinear stages of interface evolution. In addition, we detect favored tip-splitting in the shear 
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thickening case, a relevant nonlinear behavior not previously reported in the literature. We have also verified that, 
for a given mode n, smaller values of the surface tension parameter 7 lead to enhancement (suppression) of splitting 
events for 6 < (S > 0). 




FIG. 4: Snapshot of the fluid-fluid interface (t = 30 s), for the interaction of two cosine modes n = 4 and 2n = 8. Three 
values of the non-Newtonian parameter are considered <5 = (black solid curve), 8 = +0.01 (cm/s)" 2 (black dashed curve), and 
5 — —0.01 (cm/s) -2 (gray dashed curve). Other physical parameters are: a„(0) = 0.001 cm, Ro = 0.3 cm, 7 = 0.025 cm 3 /s, 
and Q — 3tv cm 2 /s. Splitting is favored for shear-thickening fluids (gray dashed curve). 



A physical mechanism that seems to be at work in shear thinning case has been proposed in references JUl [12], [13| : 
viscosity is lowered in high-shear regions, which are in front of the fingers. Hence the tip of a finger experiences less 
resistance than the part of the interface near the tip. This suggests that tip-splitting would be suppressed if the 
fluid shear-thins. Our second order mode-coupling results indicate that a similar mechanism is at work in the shear 
thickening case: the resistance would be increased at the finger tips, suggesting that tip-splitting would be increased. 



2. Action of the parameter S on side branching 

Another relevant non-Newtonian effect that can be studied at second order refers to the side branching phenomenon. 
In the framework of a mode-coupling theory, side branching requires the presence of mode 3n. If the harmonic mode 
&3n is positive and sufficiently large, it can produce interfacial lobes branching out sidewards which we interpret as 
side branching. 
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Consider the influence of the fundamental mode n, and its harmonic 2n, on the growth of mode 3n. The equation 
of motion for the cosine 3n mode is 

a 3n = A(3n) a 3 „ + i S(3n) a„a 2n , (34) 

where the side branching function jS(3n) = [T(3n, n) + T(3n, 2n)] can be easily obtained from Eq. ((33|). By analyzing 
Eq. (^) we observe that mode 3n can be spontaneously generated due to the driving term proportional to a n a2 n , 
such that it enters through the dynamics even when it is missing from the initial conditions. The existence and phase 
of mode 3n depends on the interplay of the modes n and 2n. Side branching would be favored if a 3n > 0. 




FIG. 5: Variation of the side branching function S(3n) — T(3n,n) + T(3n, 2n) with the non-Newtonian parameter 8, for modes 
n = 4,5,6. The injection rate Q = 2n cm 2 /s and 7 = 1/100 cm 3 /s. The units of 5(3n) and 5 are (cm s) _1 and (cm/s)~ 2 , 
respectively. 



To study the growth of mode 3n as the non-Newtonian parameter 5 is varied, we plot S(3n) as a function of 6 in Fig. 
5. We consider a particular R and R combination that corresponds to the onset of growth of mode 3n (i.e. obeying 
the condition A(3n) = 0) in the Newtonian limit 6 = 0. From Fig. 5 one can ve rify th at in the weak shear limit we 



consider in this work, S(3n) is negative for all values of 6. As shown in section HI B 1 , starting with a fundamental 



mode a n , the harmonic mode am is driven negative. Hence the product S(3n)a n ci2n in Eq. (B4j) is positive, driving 
a 3n > 0. 

Whether side branching actually occurs depends on the magnitude of 03^. Indeed, our discussion shows the presence 
of 03,1 > even in the case of Newtonian fluids, where side branching does not occur. Why does side branching occur 
for shear thinning fluids? We speculate its origin may lie in the cubic form of our modified Darcy's law, Eq. (|^), 
which will add a term of the form Sa^ to d 3n . Unfortunately, we have not calculated this term since we stopped 
our derivation at second order. Third order calculations for non-Newtonian Hele-Shaw flows lead to several new and 
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complex mode-coupling terms, whose analysis and interpretation are not at all obvious, and go beyond the scope of 
our current work 

Furthermore, there is a delicate interplay of mode 3n and mode 2n, described by the evolution equation ( |34| ) and 
a similar expression for the growth of mode 2n, 

a 2n = A(2n) a 2n + ^ [T(2n, n) a 2 n + S(2n) a n a 3n ], (35) 

where S(2n) = [T(2n, — n) + T(2n, 3n)] > 0. Hence, side branching via a positive a 3n will tend to drive a 2n positive 
(or at least less negative), inhibiting tip splitting but also reducing the growth rate of a 3n itself. 



C. Rectangular geometry limit 



It is interesting to study how the main interfacial features examined in the radial geometry behave if we consider 
flow of non-Newtonian fluids in rectangular Hele-Shaw cells. A mode-coupling equation describing the system can be 
obtained if we take the "rectangular geometry limit": R — > oo and Q — > oo, such that Q/(2nR) = Voo and n/R = k 
remain constant, where Woo is the flow velocity at infinity and k denotes the wave number of the disturbance. In this 
limit the interface evolution reverts to the evolution of the rectangular flow geometry, with linear growth rate given 
by 



A(fc) = \k\ [ Voo (l + Svl) - j\k\ 2 (l + 2*0] 
In addition, the only nonzero second order mode-coupling terms are 



(36) 



G N {k,k') = \k\ l-sgn(jfe k) 



(37) 



G NN (k,k') = -|fe|<4 l-sgn(fc k) 



(38) 



and 



J NN (k, k') = {3 - sgn[fc'(fc -k')]+2 sgn{k(k -k')]} 



(39) 



Note that all rectangular limit expressions (|36|)-(|39|) are time-independent, and preserve the reflection symmetries 
shown in Eq. (|30|). Based on these findings, and using Eq. (|6|) we obtain a more compact mode- coupling equation 
for non-Newtonian flow in rectangular Hele-Shaw cells 



(k = \(k)(k+ [G N (k,k') + S G NN (k,k) 

k'^kfi 

+ ^ E J NN (k,k) Cfe'4-fe'- 



(40) 



The rectangular mode-coupling Eq. (^) is useful to study the influence of weak shear effects on the development 
of the fluid-fluid interface in rectangular cells. First, we examine finger tip-splitting related issues (finger narrow- 
ing/widening). As discussed in section [II Bl, we analyze the influence of the fundamental k on the growth of its 



harmonic 2k. The equation of motion for the cosine mode 2k is 



d 2fe = A(2£:) a 2k + - T(2k,k) a 2 k , 



where 



T(2fc,fc) = ( 5[2fc 2 t- 0O (« oo - 7 fc 2 ) 21 



(41) 



(42) 



Note that the driving force term T(2k, k) vanishes in the Newtonian limit 6 = 0. This agrees with the results obtained 
in Ref. |l4j for Newtonian rectangular flow. Notice further, that for the non-Newtonian case (5 ^ 0) the harmonic 
mode 2k grows spontaneously even if missing from the initial conditions. The selected sign of a 2 k is given by the sign 
of T(2k, k) and thus is dictated by 8: if 5 > (S < 0), T(2k, k) is positive (negative), and hence a 2k > (a 2k < 0). 
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Considering the case 6 > 0, this means that the fingers become narrower as the shear thinning behavior becomes 
more pronounced. Moreover, by inspecting Eq. (^) we observe that the width of the finger tips decreases for 
increasingly larger values of the flow velocity Voq. Of course, these effects are minimized for larger values of the surface 
tension parameter 7. All these results are in agreement with recent experimental fl6|| and theoretical investigations |l7| , 
[l8t for weak shear-thinning flows in rectangular geometry, in which shear-thinning narrows the width of steady state 
fingers. 

If we take k — k* « ^Jv oa /'i r ) as the fastest growing mode, then \{2k) < so exponential growth of a 2 k is prevented. 
Unlike the case of radial flow where any mode n eventually goes unstable for sufficiently large R, the growth rates 
X(k) in the rectangular geometry are time independent. Instead of exponential growth of a 2 k, we instead expect its 
amplitude to saturate at the value T(2k, k)al/(— 2\{2k)) obtained by setting a 2 k = in Eq. (41). Alternatively, one 



could specially prepare an initial condition with an initial perturbation of wavevector k sufficiently small that 2k still 
lies within the band of unstable modes (up to about \/v 00 /'y), in which case mode 2k will be spontaneously generated 
and able to grow to a considerable magnitude. 

For completeness, we also discuss the shear thickening case: for 5 < the function T(2k, k) is negative, favoring 
finger tip-widening. Even though this shear thickening behavior has not been studied experimentally, the broadening 
of the fingers in shear thickening flow in rectangular cells has been theoretically predicted by numerical simulations [ p"8| . 
Our analytical results reinforce the correctness of such numerical predictions. 

Finally, we briefly discuss side branching behavior in rectangular Hele-Shaw cells. We consider the influence of 
modes k and 2k on the evolution of mode 3fc. In this case, the relevant equation of motion has the form 

a 3 k = A(3fc) a 3k + ^ S(3k) a k a 2k , (43) 

where, as it was in the radial flow case, the coefficient S(3k) is the sum of T(3k, k) and T(3k, 2k) 

S(3k) = S [8k 2 Voo { Voc - 7 k 2 ) ( Voo - 4 7 fc 2 )] . (44) 

The favored sign of a^k depends on the sign of the product a k a,2k and also on the sign of S(3k). If S is positive (shear 
thinning) then S(3k) is negative in the band of wavevectors 

< k < (45) 

and positive outside this band. As we saw previously in the case of mode 2fc, if the fundamental mode is taken as the 
fastest growing mode k — k* , then exponential growth of the harmonic 3k is inhibited because A(3fc) < 0. Instead it 
will saturate at a magnitude of S(3k)aka 2 k/(— 2A(3fc)). 

However, for a specially prepared initial condition perturbed at wavevector k « k* /3, the harmonics 2k and 3k will 
both be able to grow. By Eq. ( pd| ) a positive harmonic mode a 2 k > will appear spontaneously. Then, by Eq. (|43|), 
modes et^ and a 2 k will conspire to create mode a^k- Because k w k*/3 lies outside the wavevector band (|45|), the 
value of S(3k) is positive. Hence the mode a 3 k is driven positive, the sign needed to create side branches. While this 
conclusion is only tentative (due to our neglect of the complete set of third order terms) it would be interesting to 
test via numerical simulation or experiment. 



IV. CONCLUDING REMARKS 



Visually striking patterns arise when a less viscous Newtonian fluid displaces a more viscous non-Newtonian fluid 
in the confined geometry of a Hele-Shaw cell. These complex patterned structures are the result of finger tip-splitting 
and dendritic interfacial instabilities. Due to the complicated dynamics of the system, and also from the limitations 
imposed by purely linear analysis, the majority of the theoretical studies in this area of research rely heavily on 
sophisticated numerical methods. In this work, we developed an alternative theoretical approach to the problem which 
allowed us to address important nonlinear issues analytically. A key point in our derivation was the introduction of 
velocity vector potential (as opposed to scalar potential) which was demanded by the non-Newtonian fluid flow. 

We started our investigation by analyzing the role of fluid viscosity anisotropy on the development of the Saffman- 
Taylor instability in radial Hele-Shaw cells. To approach the problem analytically we considered the weak shear limit, 
and focused on the onset of nonlinear effects. In order to examine the influence of shear thinning and shear thickening 
on the shape of the emerging patterns, we derived a mode-coupling equation which is ideally suited to describe the 
weakly nonlinear interface evolution. Our analytical results show that finger tip-splitting is enhanced (diminished) in 
the case of shear thickening (thinning) fluids. We applied the rectangular geometry limit to our radial mode-coupling 
equations of motion, and also studied finger-tip and side branching behavior in rectangular cells. Our main findings 
predict finger tip-widening (narrowing) for shear thickening (thinning) fluids. 
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In neither the radial nor the rectangular geometry limit were we able to show the spontaneous appearance of side 
branching, except for some specially prepared initial conditions in rectangular geometry. We speculate that extending 
the current approach to third order could shed more light on this problem. 

In summary, our analytical mode-coupling approach detects and predicts several features of the patterns formed in 
non-Newtonian Hele-Shaw flows in both radial and rectangular geometries. Furthermore, it predicts some behaviors 
not yet investigated in the literature, such as the evolution mechanisms involving tendency towards advanced finger 
tip-splitting and reduction of side branching for shear thickening fluids. We hope the main results presented in this 
paper will prompt further theoretical and experimental work on non-Newtonian Hele-Shaw flows, especially in the 
area of early stage transients in rectangular geometry where our theory makes specific predictions that have not been 
subject to experimental or numerical check. 
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APPENDIX: SECOND ORDER MODE-COUPLING TERMS 



This appendix presents the expressions for the second order Newtonian (N), and non-Newtonian (NN) mode- 
coupling coefficients which appear in Eq. (|2£ 



and 



where 



F N (n, m) 



|n| I R 
~R I R 



1 



sgn(nm) 



7 / nm 



3m 2 



(A.l) 



F NN (n, m) 



-R 



62. 



R 



R I R 

7 2 |n| 



2a(m) \n\ 2/3(n) f(n,m) 



m n \n\ 2(\n\ + 1) 



nm 3m 2 

i? 3 (M + i) 



(A.2) 



G N (n, m) = — 



1 — sgn(nm) — - — ■ 
n 



G NN (n,m) = -rrR 2 



2a(m) \n\ 2/3(n) f(n,m) 



\n\ 2(M + 1) 



(A.3) 
(A.4) 



H NN (n, m) 



jnj R 2 
2(M + 1) R 



f{n,m) 



(A.5) 



J NN (n,m) 



2(\n\ + 1 



-R f(n,m) 



(A.6) 



f(n,m) = 



— h (n, m) + k (n, m) 
n 



(A.7) 
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